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ABSTRACT 

We test how accurately the smoothed particle hydrodynamics (SPH) numerical 
technique can follow spherically-symmetric Bondi accretion. Using the 3D SPH code 
GADGET-3, we perform simulations of gas accretion onto a central supermassive black 
hole (SMBH) of mass IO^Mq within the radial range of 0.1 — 200 pc. We carry out 
simulations without and with radiative heating by a central X-ray corona and radiative 
cooling. For an adiabatic case, the radial profiles of hydro dynamical properties match 
the Bondi solution, except near the inner and outer radius of the computational do- 
main. The deviation from the Bondi solution close to the inner radius is caused by the 
combination of numerical resolution, artificial viscosity, and our inner boundary con- 
dition. Near the outer radius 200 pc), we observe either an outflow or development 
of a non-spherical inflow unless the outer boundary conditions are very stringently 
implemented. Despite these issues related to the boundary conditions, we find that 
adiabatic Bondi accretion can be reproduced for durations of a few dynamical times 
at the Bondi radius, and for longer times if the outer radius is increased. In particular, 
the mass inflow rate at the inner boundary, which we measure, is within 3 — 4% of 
the Bondi accretion rate. With radiative heating and cooling included, the spherically 
accreting gas takes a longer time to reach a steady-state than the adiabatic Bondi 
accretion runs, and in some cases does not reach a steady-state even within several 
hundred dynamical times. We find that artificial viscosity causes excessive heating 
near the inner radius, making the thermal properties of the gas inconsistent with a 
physical solution. This overheating occurs typically only in the supersonic part of the 
flow, so that it does not affect the mass accretion rate. We see that increasing the 
X-ray luminosity produces a lower central mass inflow rate, implying that feedback 
due to radiative heating is operational in our simulations. With a sufficiently high 
X-ray luminosity, the infiowing gas is radiatively heated up, and an outflow develops. 
We conclude that the SPH simulations can capture the gas dynamics needed to study 
radiative feedback provided artiflcial viscosity alters only highly supersonic part of the 
inflow. 

Key words: accretion — galaxies: nuclei — hydrodynamics — methods: numerical 
— radiation mechanisms: general. 



1 INTRODUCTION 

Accretion of gas onto supermassive black holes (SMBHs) 
at the centers of active galaxies and the resulting 
feedback from them have strongly influ enced the for- 
mation and evolution of galaxies (e.g.. ISalpeterl 
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trends suc h as the central SMBH-host galaxy correla - 
tions (e.g., iMagorrian et all Il998l : iGebhardt et al.l I2OO0I '). 
The concordance model of structure formation based on A 
cold dark matter cosmology invokes feedback from SMBHs 
as a crucial ingredient to self-regulate galaxy and SMBH 
growth, as have been studied in numerical hydrodynamic 
simulations (e.g.. [Ciotti fc OstrikeJ [2OOII : lOi Matteo et al.l 
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Novak efal] l2010l ). as wel l as in semi-analytic models 



of galaxy format i on (e.R., Kauffmann fc Haehneltl 
Bower et al. "200^; 'Cr oton et all l2006l : I Malbon et alJ 



200c 



2007 



Somorvillo ct al. 2008). 

In galaxy formation simulations that can resolve kpc to 
lOO's of pc scales, the accretion flow onto a SMBH occurs 
in unresolved scales, and a subgrid prescription is needed to 
model SMB H accretion. The Bondi-Hoyle-Lyttlcton param- 
eterization (|Hovle fc Lvttlet on 1939; Bondi fc Hoylc 1944; 
lBondilll952l ). relating the SMBH mass accretion rate to the 
resolved larger scale properties of the ambient gas, is of- 
ten assumed in an ad hoc manner. Such an assumption has 
been made in studies of isolated systems and merg ers (e.g., 
Springel et al.. 2005: Johansson ct al. 200 9b; Lusso fc Ciottil 



2OI0I) and cosmologica l simulations (e.g..lSiiacki et al. 20071: 
Pi Matteo et al.ll2008l: [Booth fc Schav3l2009l: iDubois et all 



2OIOI ). The lack of resolution introduces uncertainties in es- 



timating the gas density and temperature near the BH and 
hence the BH accretion rate. As a numerical correction, 
the accretion rate inferred from simu lations is multiplied 
by a fact or of the order of 100 (e-g-. ISpringel et al.ll2005l: 



Siiacki et al. 2007;: iDi Matteo et al.ll2008l : iBooth fc Schavd 
20091 : ISiiacki et ai1l2009l 'l. The ISM gas density and temper- 



ature are estimated by smoothing on the scale of the com- 
putational resolution at the location of the BH. This results 
in artificially low densities and high temperatures compared 
to the case where the multiphase structure of the ISM could 
be spatially resolved on the scale of the Bondi radius. As 
a co mpensation the Bondi accretio n rate is multiplied b y 
100 ll Johansson et al.ll2009al ) or 300 (|Khalatyan et al.ll2008l ). 
iBooth fc Schavell|2009h " used a density dependent multiplica- 
tion factor. Most of these simulatio ns were performed using 
the GADGET code (|Springelll2005l ). 

GADGET is a smoothe d particle hydrodynamics (SPH , 
iGingold fc MonaghanI Il977l : ILucvI Il977l : iMonaghanI Il992l ') 
code. SPH uses a Lagrangian method for modeling fluid 
dynamics, which can handle regions of higher density with 
higher resolution, and can simulate large dynamic ranges. 
SPH is widely used in structure formation studies on cosmo- 
logical scales. The GADGET code has been demonstrated 
to pa ss several tests in the original publication (^Spri ngel 
l2005f ): standard Sod shock-tube, self-gravitational adiabatic 
collapse, isothermal collapse, interaction of a strong shock 
with a dense gas cloud, dark matter halo mass function 
and clustering, and formation of a rich galaxy cluster. 
Other studies have also used GADGET for test problems 
in various context: gravitatio nal collapse and fragmenta - 
tion of molecular cloud cores l|Arreaga-Garcia et al, 200^ 
bubble growth during hydrogen reionization I Zahn et al.l 
[2003), detection of shock waves in cosmological simula- 
tion s (IPfrommer et al.l 2006h . rotation curves of disk galax- 
ies (|Kapferer et al.ll2006l), turbulent gas motions in the in - 



tracluster medium ( Dolag et al.l l2005l: | Vazza et al] 



and baryon physics in galaxy clusters l|Borgani et al.l 



200e), 



2006f ) 



However, in spite of several studies of SMBH feedback as- 
suming Bondi accretion using GADGET, a rigorous test of 
the simple Bondi accretion, or the presentation of results 
from such a test, has been absent from the literature. 

Here, we present the results of our test simulations of 
SMBH accretion using the GADGET code. We focus on 



purely spherically symmetric accretion neglecting any mo- 
tion of the BH with respect to the ambient gas. The ques- 
tion that we want to answer is, "What is required of an SPH 
code to reproduce the Bondi problem?" . We check how well 
the code can reproduce the Bondi mass accretion solution, 
and investigate the dependence of the central mass accre- 
tion rate on parameter values, length and time scales. We 
also attempt to study one of the modes of AGN feedback: 
radiation heating. 

A rigorous estimate of BH accretion rate requires a spa- 
tial resolution of the order of a sonic radius, or at least of the 
order of the Bondi radius, which is not possible in current 
cosmological simulations. An alternative approach would be 
to perform simulations by reducing the size of the computa- 
tional volume. We follow this alternative approach, and in 
our current study we resolve the Bondi and sonic radii. Ef- 
fectively we test if it is possible to refine the subgrid model 
of SMBH accretion by measuring the BH accretion rate in 
higher-resolution simulations using the same code that is 
used for cosmological simulations. 

Similar work has been done using other numerical tech- 
niques, to which the results of our SPH- based study can 
be compared. A contemporary work by iKurosawa et al.l 
(2009) performed two and three dimensional radiation hy- 
drodynamic simulations of relatively large-scale outflows 
0.01-1 pc) fro m SM B Hs using the Eulerian code ZEUS 
Stone fc N orman 1992al lbt [Stone et al.lll993 : IClarkdll996l : 
iHaves eF^l.ii2006i ). Such a study resolving sub-pc scales can 
be used to construct a subgrid model of AGN that could 
be us ed in large-scale galaxy simulations. iKurosawa et al.l 
(|2009l ) measured the energy, momentum, and mass feedback 
efficiencies due to radiation from AGN, and found that the 
values are smaller than what is often assumed ad hoc in 
cosmological simulations. 

The foundation for the present work is based on 
the tests of the ZEUS code, which were hydrodynamic 
simulatio ns of spherically-symmetr ic Bondi accr etion onto 
SMBHs jP roga fc Bcgcl manI l2003al : iProgal l2007l ). See also 
IPark fc Rico tti (2010) for an intermediate mass BH. Once 
a code has passed the test of successfully reproducing 
the simplest case of spherical Bondi accretion, it can be 
applied to related problems in various systems. Other 
stu dies ha ve included additioiial physics to the problem; 
e.g.. |Proga fc Begelman (2003allb|) incorporated rotation and 
magnetic fields, while iMoscibrodzka fc Progal (I2OO8I . |2009| ) 
explored the effects of adiabatic index and gas temperature 
on hydro- and magneto-hydrodynamic accretion flows. 

Other author s have also u sed Bondi accretion to test 
various codes (see jEdgadl2004 for a review). In a series of 
papers, iRuffertI (|1994| ) investigated the hydrodynamics of 3- 



D classical Bondi-Hoyle accretion using a grid-based code. 
Bondi-Hoyle-Lyttleton accretion have been studied in dif- 
ferent systems, e.g., onto a protoplanetary disk using SPH 
(|Moeckel fc Throopj [2OO9I ). for binary black hole mergers 
in gaseo us environments u sing general-relativistic hydrody- 
namics l|Farris et al.ll201Gf ). in the early Universe before the 
formation of the first stars and galaxies iRicotti 2002). 

We plan to follow the same path of first testing the 
spherical Bondi accretion scenario using SPH, before ad- 
dressing the problem of AGN fee dback. One of our long t erm 
goals is to confirm the results of IKurosawa et al.1 (|2009l ) m- 
volving the feedback efficiencies. 
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This paper is organised as follows. We describe our nu- 
merical method and simulation setup in We present and 
discuss the results in ^ We summarise and conclude in 21 



2 NUMERICAL METHOD 

2.1 Revisiting the Bondi Problem 

The problem of spherically symmetric accreti on onto acen- 
tral mass was analysed in the seminal work bv lBond^ l|l952h . 
It describes a central star at rest in a cloud of gas. The gas is 
at rest at infinity, where it is parametrised by uniform den- 
sity poo and pressure Poo- The motion of the gas is steady 
and spherically symmetrical as it accretes onto the central 
star. The increase in mass of the star is ignored, so that the 
force field remains constant. The gas pressure p and density 
p are related by the polytropic equation of state, p oc p^ , 
with the adiabatic index satisfying 1^7^5/3. 

Two equations are solved for the gas motion to obtain 
the velocity v and density as a function of radius r. First, 
the equation of mass continuity. 



M = —A-KT pv = constant. 



(1) 



where Af is the mass accretion rate. Second, the Bernoulli's 
equation, which reduces to 
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where the right-hand-side represents a Newtonian gravita- 
tional potential of the central star, which is a BH for our 
case. Several types of solution are possible (Figure 2 of lBondil 
1 19521 ): the one relevant for astrophysical accretion is the so- 
called critical solution. In this solution, gas is subsonic in 
the outer parts, passes through a sonic point, and accretes 
onto the central object with a supersonic velocity. The mass 
accretion rate for such a motion is 



Me 



with 
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Here Cs — ^ ^ksT/ {nrrip) is the sound speed in the gas 
of temperature T and mean molecular weight p. Solving 
Eqs. IJ]) and ((2| gives the density and velocity of the Bondi 
solution, which we denote as ps(r) and VB{r). 

This analysis gives a characteristic length scale, the 
Bondi radius. 



Rl 



GMbh 



The latter equality in Eq. (O is for an isothermal case. These 
equations are for a purely Newtonian gravitational potential 
(Eq. As 7 ^ 5/3, the sonic radius asymptotically goes 
to zero [Rs —J- 0), i.e., there is no relevant sonic point. 

However, for a problem of BH accretion the general- 
relativistic gravitational field differs from the Newtonian 
form at very small radii. The pseudo-Newtonian Paczynsky- 
Wiita potential (which we describe in H2.21 Eq. (SJ can cap- 
ture the relativistic effects well. For the Paczynsky-Wiita 
potential as well as for the fully general-relativistic problem 
(Bcgclman 1978), the Bondi flow with 7 = 5/3 has a sonic 
point at roughly the geometrical mean of the Bondi radius 
and th e Schwarzschild radius (see also IProga fc BegelmanI 
l2003al ). 

2.2 Model Setup 

Our simulation setup consists of a spherical distribution of 
gas accreting onto a central SMBH with a mass of Mbh ~ 
10* Mq. The inner and outer radii of our computational vol- 
ume are chosen such that the Bondi and sonic radii lie well 
within our simulation domain. We choose the inner radius of 
Tin ~ 0.1 pc, which is at least an order of magnitude smaller 
than the values of Rs we explored. The outer radius is var- 
ied depending on the other model parameter values, and we 
explore a range of rout = 5 — 200 pc. 

We use the 3D Tree-PM Smoothed Particle Hydrody- 
nami cs code GADGET-3 (originally described in ' Springell 
|2005h . There are only gas particles in our simulations, be- 
cause our outer radius goes only up to 5 — 200 pc, and the 
dark matter density is much lower than the gas density in 
the central lO's - 100 pc of a galaxy. 

The central SMBH is represented by the pseudo- 



le pse i 

Newtonian potential given bv lPaczvnskv fc Wiital (|l98Ql 



GMbh 
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Here Rg is the gravitational radius of the BH. In our sim- 
ulations, Rg = 2.96 X 10^^ cm = 9.57 x 10"*^ pc; there- 
fore the Paczynsky-Wiita potential is essentially Newtonian 
within our computational domain. This is represented in the 
GADGET code by a "static potential" approach, with the 
following acceleration added to each particle, 

GMbh . 



apw 



(r- 



R.Y 



(9) 



(5) 



We also tested the effect of adding a galaxy bulge poten- 
tial in our simulation, which is described in tj3.1.2l In our 
simulations we set the gravitational softening length of gas 
to values in the range 0.005 — 0.02 pc. The minimum gas 
smoothing length is set to 0.1 of the softening lengths, which 
is 0.0005 - 0.002 pc. 



The location of the sonic point can be expressed analytically 

as 

Rs = [^^)Rb. (6) 

An important timescale is the sound crossing time from a 
distance Rb to the center (the Bondi time): 

Rb _ GMbh 

Cs Cs oo 



tB 



(7) 



2.3 Initial and Boundary Conditions 

We start with a spherical distribution of particles between 
Tin and rout, distributed according to the initial profiles of 
density pinit(r), velocity Winit(r), and temperature rinit(r). 
For most of our runs the initial profiles are taken from the 
Bondi solution, which is parametrised by the density poo and 
temperature Too at infinity. The initial conditions (ICs) are 
generated using an adiabatic index 7init, and the simulations 
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are run with, "^run- Most of oxir runs have 'yinit — 'Yrun (see 
Tables [T] and [21 . The values of different parameters we used, 
along with their justification are described in il3.f l and i )3.2l 
Any particle going out of our computational domain 
('"in < r < rout) is considered to have escaped the boundary, 
and is removed from the simulation. Particles going inside nn 
are being accreted into the inner boundary, and are counted 
in the mass inflow rate ( i|3.1.2[) . Effectively, we simulate a 
static sink of radius rin, which absorbs the accreting parti- 
cles. We tested some other outer boundary conditions that 
are discussed in il3.1.5l 



2.4 Radiative Heating and Cooling 

Radiation from the central SMBH is considered to be in 
the f o rm of a spherical X-ray emitting corona (e.g ., IProeal 
I2OO7I : IProea et al.ll2008l : [k urosawa fc Progal[2009al ). which 
irradiates the accretion flow. The X-ray luminosity, Lx, is 
a fraction fx of the Eddington luminosity, LEdd: 

A-KcGnipMBH 



Lx ~ fxLj, 



i/Edd ~ 



(10) 



where c is the speed of light, G is the gravitational constant, 
rrip is the proton mass, and Oe is the Thomson cross section 
for the electron. The local X-ray radiation flux at a distance 
r from the central source is 

Lx 



Fx = 



47rr2 



(11) 



The heating-cooling function is parametrised in terms of the 
photoionization parameter, ^, which is defined as 

47rFx Lx 



(12) 



where n — p/{fimp) is the local number density of gas. We 
use a hydrogen mass fraction of 0.76 to estimate the mean 
molecular weight fi. 

We include radiative processe s in our simu l ations using 
the heating-cooling function from | Proga et al] (|2000l ). The 
equations are originally from lBlondin l| 19941 ). who presented 
approximate analytic formulae for the heating and cooling 
rates of an X-ray irradiated optically-thin gas illuminated by 
a 10 keV bremsstrahlung spectrum. The net heating-cooling 
rate, £, is given by 



pC = n (Gcompton + Gx - it,;) [erg cm s 



(13) 



where each of the components are formulated below. The 
rate of Compton heating and cooling. 



ompton 



= 8.9 X lO-^'^C (Tx - 4r) [ergcm^s-^]. (14) 



The net rate of X-ray photoionization heating and recombi- 
nation cooling, 

T 



Gx = i.5xio-^^Ci/*r-^/^fi- — 

V Tx 



[erg cm s 



'^].(15) 



The rate of bremsstrahlung and line cooling, 

Lb,l = 3.3 X 10-27yl/2 

+ [1.7 X 10"^® exp (-1.3 X loVr) C^T 



-1/2 



+ 10-^*] 5 



r 3 -In 

[erg cm s \ 



(16) 



We adopt the optically thin version of line cooling in Eq. (|16[) 
by setting 5 = 1. In the above, Tx is the characteristic 



temperature of the bremsstrah lung radiation. We use Tx ~ 
1.16 X 10** K, corresponding to lBlondinI l|l994h 's value of 10 
keV. 

The heating-cooling rate is a function of ^, Tx and T. 
In the code, £ is computed for each active particle, and 
added to the specific internal energy (entropy in GADGET- 
3) equation of each particle using a semi-implicit method. 
In the entropy equation, the non-radiative terms are inte- 
grated in an explicit fashion using the simulation timestep, 
and then the radiative term is integrated using an implicit 
method. This integration methodology is the same as that 
of radiative cooling and photoionization heatin g in a cos- 
molo gical context in the GADGET code (e.g., iKatz et al.l 
Il996l ). 



3 RESULTS AND DISCUSSION 

3.1 Reproducing Bondi Accretion 

First, we perform a series of simulations of Bondi accre- 
tion (i.e., there is no radiative heating and cooling) as listed 
in Table [U We use Too = lO'' K and poo = W~^^ g/cm^, 
which are typical values at lO's of pc aw ay from SMBH used 
in AG N accretion simulations (see e.g., [Kurosawa fc Progal 
l2009bl . and references therein). Since Rs as y 5/3 
(Eq. |S] in ii2.1|l . we use 7init ~ 1-01 in order to have the 
Bondi and sonic radius well between rin and rout- Therefore 
the equation of state is almost isothermal, and the simula- 
tions are run with the same value of 7run — 1.01. For these 
parameters, the Bondi radius is at Rb = 3.0 pc, the theo- 
retical value of the sonic point is Rs — 1.5 pc, and the Bondi 
time is tn = 7.9 x 10^ yr. 

All the runs in Table[T]have rinit(r) = Too. In the Bondi 
IC runs, the initial condition is generated from the Bondi 
solution, i.e., the initial particles follow pinit(r) = ps(r) and 
i'init(r) = VB{r). In the uniform IC runs, we start with a 
constant initial density of poo and Uinit ~ 0. 



3.1.1 Particle Properties 

Figure [1] is a scatter plot showing the properties of particles 
vs. radius in a representative Bondi accretion simulation. 
Run 7, which has rout = 20 pc and the Bondi IC. The ra- 
dial component (vr) of the velocity v and the density profile 
follows the Bondi solution quite well, except near the inner 
and outer radii. A negative value of Vr represents inflowing 
mass, whereas positive Vr denotes an outflow. We do not 
show the non-radial velocity components (i.e., ve and v^) 
because they are typically 100 — 1000 times smaller than Vr- 
The temperature profile is almost isothermal at 10^ K {— 
Too), which is expected since we used 7iun ~ 1.01. Examin- 
ing the Mach number profile, we see that the gas is subsonic 
near rout, passes through a sonic point, and approaches rin 
with supersonic velocity (Mach — 6). The location of the 
sonic point (where the gas crosses Mach = 1) in the simu- 
lation is ~ 1.5 pc, consistent with the theoretical value of 
the Bondi solution {Rs, i)3.ip . The smoothing length (hsmi) 
of particles near rin is ~ 0.12 pc. It is much larger than the 
minimum gas smoothing length, which is set to 0.0005 pc in 
this run. This finite numerical resolution is partly respon- 
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sible for the deviations of the profiles from the ideal Bondi 
solution at small radii. 

The slight decrease of the density from ps(r) and the 
corresponding increase in smoothing length at r < 0.13 pc 
is an artifact of our inner boundary condition. There are no 
particles inside rin, as they accrete onto the sink. Therefore 
the smoothing spheres of the particles just outside of rin 
overlap with the sink, causing them to mi ss some neighbors , 
and the density is slightly underestimated. ISate et al.l (|l995l ) 
corrected the boundary conditions near the sink particle in 
SPH as a solution to this problem. Since this effect appears 
only at r < 0.13 pc in our simulations, we do not adopt any 
special boundary conditions for the sink. 

There is a group of particles flowing out of the compu- 
tational volume at rout, because of finite pressure gradient 
there. According to our outer boundary condition, there are 
no particles just outside rout, and particles just inside rout 
feel an outward pressure gradient and are pushed out. This 
explains the following features at r > 10 pc: Vr > 0, increas- 
ing Mach, decreasing p and increasing hsmi ■ We see that this 
spurious, unwanted outflow depends significantly on the gas 
temperature. Most of our simulations with a higher temper- 
ature Too = 10^ K (all in Table [T] few in Table ^ show 
this outfiow prominently, while the simulations with lower 
Too = 10^ K have much less outfiow. This is because a lower 
temperature at rout leads to a lower gas pressure, which re- 
duces the outward pressure gradient acting on the particles, 
and weakens the mass outflow at rout- 

The radial properties of all the runs in Table [1] look sim- 
ilar to Figure[T] The small scatter of the quantities vs. radius 
shows that the accretion flow is spherically symmetric, es- 
pecially in the inflowing regions. 

A cross-section slice of gas density is presented in Fig- 
ure [2] for Run 4, which has rout = 5pc and starts with the 
Bondi IC. It shows the X ~Y plane, with the cross-section 
through Z = Q. The velocity vectors of the gas are overplot- 
ted as arrows. This cross-section image shows that the flow 
is indeed nearly spherically symmetric. The velocity vec- 
tors are symmetrically pointing inward in the inner volume, 
where there is net inflow of particles. 

3.1.2 Mass Flux Evolution at the Inner Boundary 

The mass flux of accreting gas at a given time as a function 
of radius can be expressed as 

M(r) = ® p t; • da = ® p Vr dO,. (17) 

J S J 4,t: 

In particular, we compute the mass inflow rate at the inner 
boundary, Min,ri„, and consider it as the figure of merit to 
determine how well a simulation can reproduce the Bondi 
accretion. We sum up the mass of particles accreted inside 
rin within a certain time interval (At), then compute the 
inflow rate as: 

^in.n„W = ^ X! "^^^^ 

{'■fc<nn) 

The GADGET code uses individual timesteps for particles, 
so not all particles are active at once. Using the simulation 
timestep as At produced very spiky Min,ri„ over time, there- 
fore we sum over the masses of particles accreted during 128 



or 256 timesteps for the results presented in the following. 
If a simulation can track Bondi accretion perfectly, it would 
have ATin,ri„ = Mb (Eq. |3]). We check how closely and for 
how long this relation is satisfied for our runs. 

The mass inflow rate as a function of time is plotted 
in Figs. |3] and |4] for the runs in Table [1] The horizontal 
straight line in each panel indicates the corresponding Bondi 
accretion rate Mb — 4.6 x 10^^ g/s. 

Figure [3] shows the runs with the Bondi IC and a 
uniform IC (constant density and velocity), the first eight 
runs in Table [T] We see that the runs starting with the 
Bondi IC reach steady-state quickly, and can reproduce the 
Bondi accretion rate for some time. This time range in- 
creases with increasing outer radius of the computational 
volume. Runs with rout = 5, 10, 20,50 pc follow the Bondi 
rate for time durations of t ~ (0.6, 1.6,4,8) x 10'* yrs (i.e., 
0.76,2.0,5.1,10.1 is), respectively. The run with rout ~ 
10 pc is tested with 2 numerical resolutions, using particle 
numbers of A*' = 64^ and 128^. There is no significant dif- 
ference in the Min,ri„ vs. time, except that the N = 64^ 
run has a larger scatter of the mass inflow rate because of 
a lower mass resolution. The uniform IC runs start with 
a low Min.rjn at early times (lower than the corresponding 
Bondi IC run by more than an order of magnitude), have 
an increase, but cannot catch up with the Bondi rate in the 
runs with rout = 5 or 10 pc. With a larger outer boundary, 
rout = 50 pc, the uniform IC run eventually starts to follow 
the Bondi IC run after t = 3 x 10* yrs = 3.8 ts. 

The top panel of Figure U] gives the results of starting 
with the gas at rest (uinit ~ 0) and having different initial 
gas density proflles: pB{r), and a uniform density (poo). The 
mass inflow rates for the Bondi and uniform cases increase 
with time, reaches a maximum value (almost steady-state) 
of Mb , then faUs off after 4 x 10* yrs. 

The main reason for the drop in Min,ri„ after few xlO* 
yrs in all of our runs is that the particles are lost from the 
simulation volume, as described in i)3.1.1l and there are not 
enough particles between rin and rout to correctly represent 
the Bondi problem. Our simulations start with N particles 
(third columns in the tables) as the IC. Then according to 
our boundary conditions ( >i2.3|l . particles accreted into rin 
or moving outside rout are removed from the simulation. 
This causes the particle number (and hence the total mass) 
between rin and rout to continuously decrease with time. 
Despite this mass loss, several of our runs reached a steady- 
state for durations of a few xlO* yrs, i.e., a few Bondi times 
(Eq.[7|, which is enough for the Bondi problem we are study- 
ing. 

This loss of particles occurs mainly out of rout for the 
runs in Table [T] We already discussed this outflow in H3.1.1I 
that it is due to the high temperature T^x, ~ 10^ K. As 
Table [1] shows, the cumulative mass fraction accreted into 
rin is between 2 — 35% of the initial mass, by the end time 
of the simulation. The majority of mass (55 — 90%) move 
out of rout by the end of simulation. In order to remedy this 
problem, we tried other approaches of handling the outer 
boundary condition, which we discuss in i]3.1.5l 

The above results demonstrate that the simulations 
starting with the Bondi IC reach steady-state sooner than 
those starting with a uniform IC. For our next set of runs 
in Table [2j we start with the Bondi IC in order to reach the 
desired steady-state solution faster. 
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3.1.3 Hernquist Denstty & Bulge Potential 

In order to explore the dependence of our results on the ini- 
tial gas density profile, in Run f f we start with a different 
density: the Hernquist (f990) profile pH{r) and finit = 0. 
The density profile has a radial dependence of pnir) oc 
r~^(r + Obuige)"^, which is normalized to contain the same 
value of total mass within our computational volume as with 
the Bondi density profile. The corresponding mass inflow 
rate is shown in the top panel of Figure U starting with 
Afin,ri„ = f .4M_B and then a decrease with time, never reach- 
ing a steady-state. This is because pH(r) is steeper than 
pair), therefore more mass is concentrated in the inner parts 
of the simulation volume, causing faster accretion. 

The bottom panel of Figure |3] shows the effect of adding 
a galaxy bulge potential to the Bondi acc retion simulation 
setup. We consider a bulge potential for the llfernauisti (If 990l ) 
profile: 

r + abuigc 

where Mbuige = 3.4 x 10^^ Mq is the mass and flbuigc = 
700 pc is the scale length of a Milky- Way type galaxy 
bulge ijjohnston et al.lll996l ). In Run 12, we start with the 
Bondi IC, and add h to the simulation using a static- 
potential approach f i|2.2[) . Since flbuigo = 700 pc is sig- 
nificantly larger than the outer boundary rout = 20 pc, 
we are effectively adding a constant gravitational potential 
H — — GMbuigc/ibuigc to the simulation. We see that this 
increases the mass inflow rate with respect to the Bondi 
rate, which is expected because the particles fall inward and 
are accreted more readily. The value of Min,ri„ rises with 
time, but does not attain a proper steady-state, and reaches 
a maximum value of f .6 Mb before falling off. 

3.1.4 Mass Flux as a Function of Radius 

We also measure the mass flux of gas versus radius in order 
to examine how closely the flow follows the Bondi accre- 
tion rate. In our SPH simulations velocities and densities 
are computed for discrete particles, therefore we determine 
the mass flux as 

Mir) = ^ \ puvr,k. (20) 

The summation is over A'sph particles touching the spherical 
surface area at a radius r through their smoothing lengths. 
The net mass flux (Mnot) is obtained above if all are 
included in Eq. (|20[) . The inflow mass flux (Min) can be 
calculated by only counting particles with Vr < 0, and the 
outflow mass flux (Mout) by counting particles with Vr > 0. 

The resulting mass flux is plotted in Figure [5] as a func- 
tion of radius for the four runs in Table [T] with Bondi IC. 
In the inner regions (O.f < r < a few pc), all the mass is 
inflowing at the Bondi rate, i.e., Mnct ~ Min = Mb, and 
Mout = 0. Near the outer boundary, the inflow rate reduces 
and outflow increases toward rout • Here the net mass flux is 
dominated by this outflow, which is greater than the Bondi 
rate for rout ^ f pc. We flnd that the transition radii (rtr) 
from net inflow to net outflow varies with time for a given 
run. It starts at the outer boundary, and comes inward with 
the progress of time, as more and more particles flow out. 



In Run 7, all the gas between r = O.f and 8pc is inflow- 
ing at time 0.016 Myr. Beyond 8pc there is an outflow, but 
the inflow still dominates the mass flux until rtr ~ 12 pc. 
After rtr, the net mass flux is dominated by the outflow of 
the gas. The results for different outer radii are qualitatively 
similar to each other. 

3.1.5 Different Outer Boundary Conditions 

In an attempt to remedy the above mentioned outflow prob- 
lem, we tried some alternate approaches of our outer bound- 
ary condition, two of which were successful in reducing the 
outflow, as discussed below. With a simulation setup similar 
to ours, ICuadra et al.l (|2006l ) modeled the accretion of stel- 
lar winds on to Sgr A* using the GADGET code. They also 
flnd an outflow, at least 99% of their 'wind' (gas) particles 
escaping from the inner volume into the greater Galaxy. 

The numerical representation of the initial conditions 
used in our code does not guarantee that the Bernoulli's 
function (Eq. [2j LHS - RHS) is exactly zero everywhere. 
Therefore we check whether there are any unbound particles 
in the initial conditions, and if so what are their effects. 
We found that the fraction of initially unbound particles 
depends on the initial scatter of the radial density proflle. 
But these unbound particles have no significant effect on the 
total outflowing mass fraction. We performed a few more 
tests starting with an initial condition with < 5% scatter 
(at a fixed-r) in the Bondi density profile, and one with 
a perfectly uniform IC. These runs also have the outflow, 
with the same mass fraction flowing out of rout as the larger 
initial-scatter runs. This verifles that the outflow is being 
caused primarily by our outer boundary condition, i.e., the 
outward pressure gradient on particles near rout because of 
missing neighbors. 

We note that others also found that setting up the de- 
sired boundary con ditions can be very challenging in SPH 
(e.g.. lHeran"tHf994l ). The main reason for this is that the 
quantities, of which conditions are to be set, are determined 
by smoothing over neighboring particles, which could be lo- 
cated anywhere with respect to the boundaries. For our case 
it is the difficulty of implementing the Bondi assumption of 
an inflnite and spherical reservoir of gas at inflnity. In an Eu- 
lerian code like ZEUS the gas density at the outer boundary 
can be easily held fixed at a constant value, therefore an 
outer radius as small as 1.2R b is enough to adequately re- 
produce Bondi accretion (see lProga fc Begelman|[2003a^ . 

A similar outflow issue is dealt with in MHD SPH simu- 
lations of molecular cloud core collapse by placing the cloud 
within a uniform, low-density box of surrounding material 
in pressure equilibrium, where the outer box size is twice the 
cloud radius in each direction (|Price fc Bate|[2007l ). 

We are successful in reducing the outfiow to a negligible 
fraction, still maintaining spherical-symmetry of the accret- 
ing gas, by altering some property of every particle which 
has a net outward radial velocity {vr > 0) in a wide outer 
shell: between 0.4rout and rout, in a run with rout = 20 pc. 
Below we describe the runs with two of the stringent condi- 
tions which can successfully reduce the outflow. Runs 7a and 
7b have everything similar to Run 7, but with two different 
implementations of the outer boundary condition. 

(f) In Run 7a, whenever a particle between [0.4rout — 
rout] has Vr > 0, we change its velocity to the Bondi veloc- 
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ity, i.e., we set Vr — vb, keeping everytliing else (r, 9, cj>, etc.) 
same. Since vb is negative, the outward flowing particle is 
forced to move inward, on imposing this condition. As listed 
in Table 1, such a boundary condition can reduce the out- 
flowing mass significantly: 65% mass is accreted into rin and 
4% flows out of Tout, within lOOts. 

(2) Run 7b has a layer of particles at rout which do not 
feel the net outward pressure gradient because of missing 
neighbors just outside rout. We alter the (VP/p)i term di- 
rectly in the particle acceleration dvi/dt equation. We set 
a zero pressure gradient, {VP/p)i — 0, for all particles be- 
tween [O.Qrout — rout] and for those which have t^r > in 
between [0.4rout — rout]. This run has a further reduction in 
the outflowing mass: 93% accretes in, and only 1% moves 
out, within 125ts. 

Figure [6] shows the mass inflow rate as a function of 
time of Runs 7a and 7b. We see some initial unsteadiness 
in both the runs manifesting as spikes in Min,ri„. After 0.08 
Myr, the mass inflow becomes nearly-steady (~ AIb) for 
~ 12t_g in Run 7a and ~ 25ts in Run 7b. Min.r^^ decreases 
after that; now this reduction is because a dominant mass 
fraction is being accreted into rin. 

Such an outer boundary condition (as in Run 7b) can 
be adopted with any of the runs presented in this paper 
to obtain a longer duration for which steady-state Bondi 
accretion continues onto rin. All the gas properties remain 
the same as presented here, for that longer duration. 

Table [3] shows some relevant time durations for the 
Bondi IC simulation runs, with increasing rout. The ra- 
tio of the total initial mass to the Bondi accretion rate, 
Mtot,ic/ Mb , gives an upper limit to the time duration for 
which the simulation gives reliable results. Mtotjc / Mb is 
the time within which all mass would be accreted in, assum- 
ing accretion at the Bondi rate. In practice the significant 
time duration is rather a fraction of Mtot.ic/ Mb, as we see 
in our runs, because there has to be some particles remaining 
in the simulation to reproduce the Bondi problem. 

Summarising the results, our simulations can reproduce 
adiabatic Bondi accretion within a limited time duration, 
ranging between one to ten Bondi time, which is long enough 
to investigate the Bondi problem. The radial profiles of gas 
density, velocity, temperature, and Mach number match the 
Bondi solutions quite well, except near the inner and outer 
radii. 



3.2 Spherical Accretion with Radiative Heating 
and Cooling 

In the next series of simulations we include radiative heat- 
ing and cooling. These runs are relevant to the studies of 
feedback due to radiative heating in the context of SMBH 
accretion. Table [2] gives a list of the runs, where we follow 
the prescriptions in i]2.4l and our spherical accretion sim- 
ulation setup. All these runs are done with 7run = 5/3. 
Since Rs — Opc for 7run ~ 5/3, the initial condition is 
generated using the Bondi solution with 7init = 1.4 in or- 
der to have the sonic point within current computational 
volume. We varied the X-ray luminosity in the range of 
Lx /^Edd = 5 X 10~^ — 0.5. We also vary several other simu- 
lation parameters in our heating-cooling runs (rout, N, 7init, 
Too, Poo, Tinit), as dlscussed in the subsequent sections. 

The radiative equilibrium condition can be obtained by 



solving for C — in Eq. (I13|l . This gives a relation be- 
tween the equilibrium temperature. Trad, and the photoion- 
ization parameter. The Compton equilibrium temperature 
is Tc = Tx/4 = 3 X 10^ K, using Tx from gg^] This is com- 
parable to the equilibrium Comp ton temperature of 2 x 10^ 
K found in|s azonov et al.l (|200J). At high photoionization, 
^ > 10^, Compton processes dominate, resulting in a con- 
stant temperature Trad ~ Tc for all ^. At ^ < 10^, other 
radiative cooling components start dominating and Trad de- 
creases below Tc- 



3.2.1 A Representative Case 

Run 23 is a representative simulation with radiative pro- 
cesses included, and the particle properties as a function of 
radius is plotted in Figure [T] The velocity and density pro- 
files are qualitatively similar to that in the Bondi accretion 
Run 7 in Figure [1] the gas fall inward with a larger velocity 
at smaller radii and its density rises. 

The temperature profile in Figure [7] shows signatures of 
heating and cooling operating at different radii, in contrast 
to the isothermal temperature profile of the Bondi problem. 
Within r ~ 0.7 — 180 pc, temperature is between T 10* — 
10^ K and slowly increases inward. There is a sharp increase 
of temperature toward rin, and a slight heating near rout- 
The sonic point is close to the outer boundary at rout = 200 
pc, and the gas in most of the volume is supersonic. The 
Mach number increases inward up to ~ 30 at r ~ 0.7 pc, but 
decreases at smaller r, as opposed to the monotonic increase 
in the Bondi accretion run. This is because of the steep 
heating near rin in Run 23, which causes a greater increase 
in the sound speed than the increase in the dynamical gas 
velocity, resulting in a decrease of the ratio \v\/cs. 

Supersonic gas (with 7run = 5/3) should only feel the 
gravity of the central SMBH effectively in this run, because 
the gas pressure would be negligible compared to the central 
gravity; therefore the gas should essentially be in the state 
of a free-fall. We compute the theoretical free-fall scaling of 
some of the gas properties, and compare those with the par- 
ticle properties in our numerical simulations. The free-fall 
velocity for radial infall (in a purely Newtonian 1/r poten- 
tial) is 



I2GMbh 
Vff = -\l r. 



(21) 



Then taking the equation of continuity, assuming steady- 
state, only-radial dependence of density, and using the above 
vff, we obtain the free-fall density scaling, 

Pll = (l-Y"\ (22) 
po \roJ 

Here, po and vq are scaling constants, and the density is 
equal to po at r = ro. Using Vff and pff we obtain a free- 
fall adiabatic radial profile for the temperature by solving 
the equation of internal energy with only the adiabatic term 
and 7 = 5/3, 

Tff.a « r-\ (23) 

We also solve the general internal energy equation including 
both adiabatic and radiative terms, using «// and pff, and 
obtain the temperature solution Tff^ar- The equations and 
details of the method are given in Appendix lAl Mach//,ar 
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is obtained by calculating the sound speed from Tfj\ar, and 
taking its ratio with the free-fall velocity. 

The free-fall scalings Vff, pff, Tff^ar and Mach//_a,. 
are overplotted as the dashed curves along with the particle 
properties in Figure [T] The values of the scaling parameters 
are manually set to po = 3 x 10"^'^ g/cm'^ and ro = 65 pc, 
in order to get a good agreement. We see that the particle 
properties in our simulation match the theoretical free-fall 
solutions well in the intermediate radii, as predicted for su- 
personic gas. There is a deviation near the inner radius, 
where the gas velocity in our run is lower than Vff, for 
which numerical resolution is partly responsible. Particles 
are heated up more than Tff^ar at r < 0.7 pc, and have a 
steep rise in temperature reaching T — few x 10^ K at rin, 
compared with the expected value of few x 10^ K. This is 
reflected as a decrease in Mach number at r < 0.7 pc, that 
drops to ~ 2 at rin, which is much smaller compared to the 
expected value of 30. The region near the outer radius also 
shows a deviation from the free-fall scalings. This is caused 
by particles flowing outward at rout. 

In part the heating occurring near the inner boundary 
is caused by adiabatic compression of the gas flowing inward 
into smaller volumes. The free-fall temperature relation ex- 
pected from adiabatic processes only (Tff^a) is plotted as 
the dotted curve in Figure [T] Comparing the temperature 
between the simulation and two predicted radial profiles, we 
see that the slope of T vs. r of the simulation particles near 
the inner radius is larger than both the predictions Tff^ar 
and Tff^a- Therefore this inner heating cannot be fully ex- 
plained by adiabatic and/or radiative processes. On further 
investigation (see below for more details) we find that this 
extra heating of the gas is caused by the artificial viscos- 
ity (AV) in the GADGET code. In Run 23, the AV heating 
dominates at r < 0.7 pc. The extra thermal energy comes at 
the expense of kinetic energy. The velocity of the particles 
is smaller than Vff at the same r, as the top-left panel of 
Figure [7] shows. 

The AV heating occurs within a radius much smaller 
than the sonic or critical radius, which is close to the outer 
boundary rout ~ 200 pc. Therefore the gas being heated 
up by AV is already supersonic, and remains so as it falls 
into rin. The Mach number reduces, however the gas does 
not become subsonic near rin, as the bottom-right panel in 
Figure [7] shows. So the gas temperature is affected by AV, 
but the central mass inflow rate is not affected because the 
gas is still supersonic. 

We perform 1-D simulation with the ZEUS code, us- 
ing the same parameters as Run 23, for comparing with 
the GADGET run. The particles properties from ZEUS are 
overplotted in Figure [7] as the solid curve in the first three 
panels. We see that the ZEUS results follow the free-fall 
scalings well up to the inner radius, hzeus is essentially in- 
distinguishable from Vff. pzBus is at a constant offset from 
Pff because of different normalization. Tzbus goes to some- 
what higher value than Tff^ar, reaching Tzbus ~ 10® K at 
rin; however not as high as the temperature in the GADGET 
run (r = few x lO'^). The difference between GADGET and 
ZEUS results of the temperature and velocity near rin fur- 
ther confirms that the inner heating is a numerical artifact. 
It is due to AV, which is only present in the GADGET run. 

The mass inflow rate at the inner radius of Run 23 is 
plotted in Figure [8] (bottom-right panel). A/in,ri„ increases 



with time reaching a value 4 Mb att — 2 Myr, and continues 
to rise after that. It takes a longer time to reach a steady- 
state, compared to the Bondi accretion runs. However, Run 
23 reaches almost steady-state at later times. 

3.2.2 Dependence on Model Parameters 

We explored the effects of varying some model parameters: 
rout. A'', 7init, Too, Poo, and Tinit, in our spherical accretion 
simulations with radiative heating and cooling. The mass 
inflow rate vs. time for the first four runs in Table [2] with 
Tinit = Too = 10^ K is shown in Figure[51 Runs 13 and 14 are 
shown in the top panel, for which poo = 10~^^ g/cm'^ and the 
Bondi accretion rate is Mb{j = 5/3,poo,roo) = 4.9 x 10^* 
g/s. This panel gives results with different outer radii, and 
we find Min,ri„ (rout = 50 pc) < Mn,ri„ (rout = 20 pc) up 
to t — 3 X 10* yrs, the time until which the 20 pc run has 
enough particles within the simulation volume to track the 
Bondi problem. 

The evolution of A/in,ri„ for these four runs using Tinit = 
10^ K are qualitatively similar. The runs start with a high 
inflow rate of 1.4Ms (Run 14) and 2Mb (Runs 13, 15, 16), 
which is the maximum value they reach, and have a decline 
in A/in,ri„ with time. Runs 13 and 14 never reach a steady- 
state. These runs have a very large mass fraction moving out 
of the computational volume, the largest among all the runs 
we have done. As given in Table [21 98 — 99% of the total 
initial mass moves out of rout by the end of the simulation, 
while only 0.1 — 0.5% of the total mass is accreted into rin. 
This is because the high initial temperature Tnit ~ 10^ K 
heats up gas, which tends to expand and increase the outflow 
at rout. The mass outflow at the outer boundary is almost 
solely controlled by the gas temperature, because changing 
other parameters has almost no effect on it ( i|3.2.4[) . 

The mass inflow rate as a function of time for seven runs 
in Table [21 which has Tnit ~ T^ad (the radiative equilibrium 
temperature), is shown in Figure [8l The time dependence 
of A/in, of these runs is similar to each other; it increases 
with time, reaches a maximum value, then falls off after a 
certain time. Runs with rout ~ 20 pc have an increase in 
Mn,ri„ up to ~ 10^ yrs, and for runs with rout = 200 pc, 
Afin,ri„ continues to increase even 2 Myr into the run. 

The adiabatic index of the IC is changed from 7init = 1.4 
to 1.6667 from Run 17 to 18, and in the later runs 7init = 
1.6667 is used. This change in 7init has a relatively small 
impact on the mass inflow rate vs. time (as can be seen 
from the bottom-left panel of Figure [8]), with a larger 7init 
producing slightly lower Min,ri„ ■ 

Three runs (17, 18, 21) are done with Too = 10^ K, 
which is the same Too as our spherical Bondi runs, for which 
the theoretical Bondi accretion rate is Mb = 4.9 x 10^* 
g/s. In these simulations, the mass inflow rate is orders 
of magnitude higher than Mb, and reaches a maximum 
value of 80Ms in Runs 17 and 18, and 2000iW"B in Run 
21. The reason is that in these runs the gas temperature 
is calculated self-consistently including adiabatic and radia- 
tive heating-cooling mechanisms, which reaches a value of 
T(rout) - 10' K near the outer boundary (shown in Fig- 
ure [To] and discussed in i]3.2.3p . This outer temperature 
T(rout) ~ 10' K is lower than Too = 10^ K, and implies 
a larger mass accretion rate for the corresponding Bondi 
solution that the simulation is relaxing to. 
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Therefore we do the remaining simulations with Tod ~ 
10'^ K in Table [1 Runs 19 and 20 have = 10"^^ g/cm^, 
where the Bondi accretion rate is Mb ~ 4.9 x 10^^ g/s. As 
the top-left panel of Figure[8]indicates, the runs reach almost 
steady-state when considering time intervals of few x 10* yrs, 
having an inflow rate between 1 and 2 Mb- A larger outer 
radius rout = 200 pc produces an increasing A^in,ri„ for 0.7 
Myr (shown in top-right panel) , reaching a maximum 6Mb ■ 

Runs 22 and 23 have the same simulation parameters, 
but a different particle numbers oi N = 128^ and 256^. 
With a lower resolution, Run 22 is continued until the Bondi 
problem is represented correctly, or, as long as the mass 
inflow rate increases ( i)3.1.2p . reaching a maximum Min,ri„ ~ 
5Mb after 3 Myr. 

Comparing the T vs. r plots of Runs 22 and 23, we de- 
termine the dependence of artificial heating on resolution. 
We find that the temperature up to which the gas is ar- 
tificially heated (T ~ 10* K at rin) remains the same on 
changing resolution. However, the radius within which the 
extra heating occurs increases with lower resolution; in Run 
22 gas within 1 pc is heated, which is ~ 0.6 pc in Run 23. 
This implies that a greater volume of gas is being affected 
by AV heating as resolution degrades. 

Reducing the initial temperatures to Tinit ~ Tlad 
brought a noteworthy change in our results: a larger mass 
fraction is accreted inside the inner radius rather than flow- 
ing outside the outer boundary. In the relevant runs 90—99% 
of the total initial mass accretes into rin, and only 1 — 5% 
of the total mass moves out of rout- With radiative heat- 
ing and cooling incorporated, the gas temperature near the 
outer radius is between T(rout) ~ 10** — 10^ K, which re- 
duces the gas pressure w.r.t. our previous runs which had 
T(rout) = Too = lO'^ K. This reduced outward pressure 
gradient lowers the unwanted mass outflow at rout, and in- 
creases the desired inflow. 

3. 2., 3 The Thermal History of Particles 

In Figure [TD] we plot the time evolution of particles in the 
T— ^ and T — r planes of Run 18. Note that in the T— ^ plane 
the right-most region (highest-^) is near the inner radius, 
the middle portion is the middle volume, and the left region 
is near the outer boundary. Initially the gas is in radiative 
equilibrium, therefore the particles lie on the Trad— C curve at 
t = 0, where the inner particles has a temperature T ~ 10^ 
K. The next epoch plotted is t = 0.004 Myr, and after that 
in intervals of 0.04 Myr. 

In the T — r plane, with the progress of time up to the 
epoch t = 0.12 Myr, more and more particles at r < 10 pc 
are cooled to T = lO** — 10^ K by mostly radiative cooling. 
There is a non-radiative heating near rin causing particle 
temperature to rise to T = 10^ — 10* K. The radius inside 
which this heating occurs decreases with time; particles in- 
side r — 0.6 pc are heated at t = 0.04 Myr, which reduces to 
heating inside r = 0.15pc at t = 0.12 Myr. We also see some 
heating at r > 10 pc increasing in prominence with time. 
This is radiative heating because of the decreasing density 
near the outer boundary caused by the outflow. There is a 
cooling very close to r = rout = 20 pc caused by adiabatic 
expansion, more prominent in the later epochs. 

The above features near rout manifest themselves in the 
T — ^ plane as the portions below the equilibrium curve, 



stretched out in ^. Particles near the outer computational 
volume evolve along the T^ad — C curve from lower left to 
upper right as long as radiative processes dominate. Then 
they move out of the curve and stay at lower T, when ra- 
diative and adiabatic terms become comparable, and finally 
adiabatic cooling dominates at rout. 

By the last epoch shown in Figure [10] t — 0.16 Myr, the 
simulation has started to deviate from the Bondi solution, 
because the mass inflow rate of Run 18 (Figure [S] bottom- 
left panel) increases only up to 0.13 Myr. The density of the 
whole volume at t = 0.16 Myr is lower than the Bondi profile 
by an order of magnitude. Therefore the photoionization 
parameter is higher, and the gas is heated up more than the 
temperature in the previous epochs uniformly at all radii. 

Figure [TT] shows the temperature vs. photoionization 
parameter of four runs in Table [2] which has Tinit = Too = 
10^ K, along with the T — ^ equilibrium curve marked as 
the solid line. The gas in these runs is highly photoionized 
(10^ < ? < 10"). The T vs. ^ of nine other runs (having 
Tinit ~ Trad) in Table[2]are shown in Figure [T^ where the gas 
has a lower photoionization (10~^ < ? < 10^). Each run has 
a subset of particles falling on the T — ^ equilibrium curve, 
implying that those are dominated by radiative processes. 
The increase in T over T^ad at large-^ (upper-right in the 
T — ^ plane) is due to non-radiative heating dominating near 
rin. The drop in T below T^ad and the scatter at the left in 
Figure 1111 and that at intermediate-^ (lower-right in the 
T — ^ plane) in Figure [12] is caused by adiabatic cooling due 
to the outflow at rout. 

The non-radiative heating toward rin seen in Figs. llOllllI 
and [12] is caused by a combination of adiabatic compression 
and numerical AV near the inner radius. 

3.2.4 Dependence on X-ray Luminosity 

We explore the dependence of gas accretion properties on 
the X-ray luminosity of the corona by performing runs with 
different Lx- Runs 15 and 16 are done with Lx /L^dd ~ 0.5 
and 5 x 10~*, whose mass inflow rates are shown in the 
bottom panel of Figure [9] These runs have poo = 10~^^ 
g/cm"^, corresponding to a Bondi rate of Mb = 4.9 x 10^* 
g/s. The run with the lower Lx produces a higher inflow 
rate. Between t = (1.5 — 2.5) x 10* yrs, when the runs 
are close to steady-state, Min.ri„ is 1 and 0.8 times Mb for 
Lx/L-Edd = 5 X 10"* and 0.5, respectively. A higher X-ray 
luminosity corona heats up the gas more, causes the gas 
to expand, hence increases the outflow rate at rout and de- 
creases the inflow rate at rin. But the reduction in Min,ri„ is 
smaller only by a factor of 0.8 when Lx is increased by 10^ 
times, in this particular case. This is because of the very low 
gas density of these runs, which causes the gas to be highly 
photoionized with these Lx (Figure llip . The average pho- 
toionization parameter decreases from $,{Lx = 0.5) ^ 10^^ 
to £,{Lx — 5 X 10~*) ~ 10*, which is still high, and the gas 
remains around the 'hot-branch' of the T — radiative equi- 
librium curve. The low-density gas is almost equally heated 
up with these values of Lx- Consequently, we also see that 
decreasing Lx/L^dd from 0.5 to 5 x 10~* has no effect on 
the outflow (~ 0.98) and inflow (~ 0.005) mass fractions in 
these simulations. 

We investigate more varying-Lx cases by restarting 
Run 23 (with Lx/L^dd = 5 x 10"*) at t = 1.4 Myr, 
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and performing subsequent Runs 24 — 28 with: Lx/L^dd ~ 
5 X 10-^5 X 10-^ 0.01, 0.02, 0.05. With larger Lx, parti- 
cles have higher photoionization parameter, and shifts to 
the right of the T — ^ equilibrium curve. This can be seen by 
comparing the T — ^ planes of the six runs in Figure [T2] in 
lower two panels; particles in Run 24 (lowest Lx) occupies 
the left-most area, in Runs 23 — 28 they progressively shift 
right as Lx increases. Also with rising Lx, particles are in 
radiative equilibrium up to a higher temperature: T ~ 10® 
K in Run 25, up to T ~ Tc = 3 x lO'^ K (the Compton equi- 
librium temperature) in Run 28, because radiative heating 
dominates with such X-ray luminosity. 

The mass inflow rate at the inner radius of Runs 23 — 
28 are shown in Figure 1131 We see that Min,ri„ drops at a 
varying rate as Lx is increased, because the gas is heated 
up and expands. There is only a slight difference in the mass 
inflow rate by changing Lx from 5 x 10~^ to 5 x 10~^. The 
two lowest Lx runs (24 and 23) have the same Min,ri„, and 
Lx/L-Edd = 5 x lO^'' (Run 25) produces 0.95 times lower 
Ai/in,ri„. The central inflow rate becomes 0.8 times lower with 
Lx '/LEdd = 0.01 in Run 26. 

Further increases of the value of Lx by a factor of 2 or 
more show signiflcant changes in the gas dynamical behav- 
ior. The mass inflow rate decreases by 3 — 4 orders of mag- 
nitude, and an outflow develops because the gas is heated 
up by the higher X-ray luminosity. We see that the transi- 
tion from net inflow to net outflow occurs with a value of 
Lx/L-Edd between 0.01 and 0.02, for poo = 10~^^ g/cm'' and 
Too = lO'^ K. The mass inflow rate in Run 27 rises again af- 
ter 2.2 Myr, because the outflow is not strong enough with 
Lx /L^dd = 0.02, and the inflow resumes later. There is how- 
ever a strong outflow with Lx/L^dd ~ 0.05 which expels a 
substantial fraction of gas out of the computational volume, 
and Min,ri„ never rises within 5 Myr. 

We find that, in Runs 26 and 27, the gas cools and 
gets denser non-spherically leading to fragmentation after a 
time of ~ 1.6 Myr. Then the mass inflow rate in these runs 
becomes very noisy (Figure [T3}. This fragmentation appears 
to be related to a growth of thermal instability. However, 
numerical effects might also be of some relevance which we 
will investigate in the future. 

The temperature profile as a function of radius for 
the three lower-Lx runs are shown in Figure 1141 With 
Lx/L^dd = 5 X 10^^ and 5 x 10^*, particles in most of 
the volume have a low temperature T < 10^ K. There is a 
significant heating at r < 0.2 pc in Run 24, and at r < 0.6 pc 
in Run 23, increasing temperatures to T ~ lO'^ — 10* K. With 
Lx/LEdd = 5 X 10""^, particles at all radii are heated with 
an almost constant slope. As Lx increases, ^ rises and the 
particles reach a higher radiative equilibrium temperature, 
so the whole T — r curve shifts upward in T, visible in the 
panels from left to right in Figure [TJ] This higher T^ad and 
the inner heating combined together has an indistinguish- 
able slope in the T — r plane of Run 25. 

Indications of AV heating toward rin can be seen in 
Figure 1141 The predicted temperature profiles Tff^ar and 
Tff,a ( ^3.2.ip are plotted as the solid and dashed curves in 
each panel. The simulation particles near the inner radius 
are heated above both the predictions. It is more prominent 
in Runs 24 and 23 (left and middle panels), but all the three 
runs show it. The radius inside which the AV heating occurs 
depends on Lx- Gas inside r < 0.2 pc is artificially heated 



with Lx /LEdd = 5 X 10 ^ , and the gas inside r ^ 0.5 pc for 
Lx/LEdd ^ 5 X 10--*. 

The gas properties as a function of radius in Run 26 
{Lx/LEdd ~ 0.01) is plotted in Figure [TSl where the gas is 
radiated by 20 times higher X-ray luminosity as compared 
to Run 23 {Lx /^Edd = 5x lO"", FigureQ. The radial veloc- 
ity in Run 26 is slower than free-fall, but the density profile 
is qualitatively similar to that in Run 23. The temperature 
profile shows radiative heating by a higher X-ray luminos- 
ity, up to the inner boundary. The Mach number reduces, 
because the gas becomes hotter increasing the sound speed, 
and the gas is mildly super-sonic (Mach ~ 1.5 — 2) as it 
accretes onto rin. 

Figure [TS] shows the radial gas properties of Run 28 
{Lx/LEdd = 0.05), where the gas is violently heated and 
expands. There is a well- formed outflow between r ~ 3 — 
25 pc, where all the particles have Vr > 0. This is to be 
distinguished from the outflow near rout owing to the outer 
boundary condition. The relevant outflow in Run 28 occurs 
at intermediate r values, and is surrounded by a region of 
inflowing gas between r ~ 25 — 150 pc. So it is clearly the 
previously free-fall inflowing gas in Run 23 that has reacted 
to the increased Lx, has slowed down, and is eventually 
outflowing at the time plotted (1.5 Myr). This outflow is 
caused by the X-ray heating of the gas due to the high Lx 
value in Run 28. 

The intensity of the heating can also be seen from the 
scarcity of particles at r < 0.7 pc in Figure [16] The gas is 
heated and expands, and the rate of inflow from the outer 
volume has reduced. This causes a decrease in density near 
Tin, about 50 times lower than Run 26. The gas becomes sub- 
sonic almost throughout the volume, except Mach 1 for 
few accreting particles next to rin. The temperature profile 
shows radiative heating at all radii. The small temperature 
spike at r ~ 25 pc corresponds to the heating caused by the 
transition from outflowing gas inside to inflowing gas outside 
25 pc. 

Summarising the effect of X-ray luminosity of the cen- 
tral corona, we find that a higher Lx produces a lower mass 
inflow rate at the inner boundary. This is because of feed- 
back due to radiative heating occurring in our simulations. 
With high enough value of Lx, the gas thermally expands. 

We perform some 1-D and 2-D simulations with the 
ZEUS code for comparison. These calculations conflrm some 
of the results we obtained in our GADGET runs. In the 
ZEUS counterpart of Run 23, the mass inflow rate at the 
inner boundary settles at the value of 5 x 10^^ g/s after a 
time of 8 Myr. This time duration is 4 times longer than the 
duration for which the GADGET Run 23 is evolved. To con- 
tinue the GADGET run for a comparable time, significant 
computational resources would be required. ZEUS counter- 
parts of simulation Runs 26 and 27 also show central mass 
inflow rates consistent with our GADGET results. 



3.2.5 Controlling Artificial Viscosity Heating 



The standard AV 
(jMonaghan fc Gingol 



pres cription used in SPH 

^ dl 19831) is i mplemented in the 

GADGET code fsee ISpringell |2005| . Eqs. 9-17). The 
velocity of particles is reduced by a viscous force, which 
generates extra entropy. This transforms kinetic energy of 
gas motion irreversibly into heat. The parametrization of 
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AV in GADGET is taken from iMonaghanl ([l993). All our 
simulations are done with a viscosity parameter a = 0.8 
(values between a ~ 0.5 — 1.0 are typical). 

This AV heating is a numerical artifact in some of 
our simulations, which increases the gas temperature higher 
than adiabatic and radiative predictions. It could be present 
in any other studies using SPH, and without proper testing it 
could lead to wrong conclusions about heating mechanisms. 
In our simulations, AV does not affect the mass inflow rate 
at the inner boundary since only the highly supersonic parts 
of the inflow are heated. But one must be careful if using the 
properties of this flow to investigate gas behavior on scales 
smaller than the inner radius (0.1 pc), and infer some impli- 
cations for processes like star formation and the formation 
of an accretion disk. 

In certain simulations it might be possible to control 
the magnitude of AV, or try an alternate AV scheme to 
minimize such undesired dissipation. In general, however, 
the adverse effect of AV is unknown in an SPH simulation, 
making it h ard to disentan gle its effects from other known 
mechanisms. iThacker et aP (2000) found that AV is the sin- 
gle most important factor distinguishing the results when 
they analysed twelve different implementations of SPH in 
cosmological simulations. 

Artificial viscosity is needed in SPH to prevent inter- 
penetration of SPH particles and thus a multi-valued flow, 
to allow shocks to form, and to damp post shock oscil- 
lations. However, AV might cause unphysical dissipation 
away from shocks and unwanted heating. In order to re- 
duce these problems, various modifications to the stan- 
da rd AV r e cipe i n SPH have been proposed: the switch 
bv iBalsaral (|l99%_ using time-varying visc osity coefficient 
for each particle ( Morris fc Monaghanll997l ). restricting the 
von Neumann-Richtmeyer AV to only non-compressed gas 
and calculating the bulk AV from a n integrated mea n ve- 
locity in a particle's neighbourhood (|Selhammaill 19971 ) . us- 
ing the total time derivative of the velocity divergence as a 
shock indicator and adapting individual viscosities for par- 
ticles ijCullen fc Dehnenll2010l V Other dissipatio n-inhibiting 
modified A V prescriptions have been reviewed in IMonaghanl 
(|2005l ) andd osswod (|2009l . and references therein). 

A few tests of modified AV schemes in certain astro- 
physical scenarios show that the proposed improved AV 
recipe s are not always effective. iCartwright fc Stamatello^ 
l|2010h tested the Balsara switch and time dependent viscos- 
ity for their effectiveness in disabling AV in non-convergent 
Keplerian shear fiow, in SPH simulations of accretion disks. 
They found that neither of the mechanisms work effectively, 
as th ose methods leave AV active in areas of smooth shearing 
fiow. lJunk et all (|2010l ') showed that the Kelvin-Helmholtz 
instability is increasingly suppressed for shear flows with ris- 
ing density contrasts, where using the Balsara switch does 
not solve the problem, indicating that other changes to the 
SPH formalism are required in order to correctly model 
shearing layers of different densities. 

In a future study we will test if some alternate AV pre- 
scription can reduce the excessive inner heating occurring in 
our spherical accretion simulations. 



4 SUMMARY AND CONCLUSION 

Aiming to answer the question how well can the SPH tech- 
nique follow spherical Bondi accretion, we perform numer- 
ical simulations of the Bondi problem using the 3D SPH 
code GADGET. Our simulation setup comprise of a spher- 
ical distribution of gas, between the length scale 0.1 — 200 
pc, accreting onto a central SMBH. Represented by discrete 
particles in SPH, the gas undergoes hydrodynamic inter- 
actions, and experiences gravitational potential of the BH, 
implemented by the Paczynsky-Wiita 'static' potential. The 
values of model parameters we explore are selected in or- 
der to encompass the Bondi and sonic radii well within our 
computational volume. 

Our spherical accretion simulations (Table [1] 7run = 
1.01) can reproduce Bondi accretion within limited temporal 
and spatial domains. The main results are summarised below 
(points 1 — 4): 

(1) The radial profiles of the velocity and density of 
the simulation particles match the corresponding Bondi so- 
lutions quite well, except near the inner and outer radii. 
The temperature profile is almost isothermal at 10^ K {— 
Too). The gas is subsonic near rout, crosses a sonic point 
at ~ 1.5 pc (consistent with the predicted Rs value of the 
Bondi solution), and approaches rin supersonically. 

The deviations of the profiles from the ideal Bondi so- 
lution at small radii is caused by three effects: numerical 
resolution (smoothing length of the SPH particles), no spe- 
cial c orrections to the inner boundary conditions (Bate et al.l 
Il995l ). and artificial viscosity in the SPH code we are us- 
ing. The deviations at the outer boundary is caused by the 
outflow of particles, occurring in several of our simulations, 
because of the outward pressure gradient at rout. 

(2) Runs starting with the Bondi IC reach steady-state 
quickly (earlier than starting with a uniform IC). The mass 
inflow rate at the inner boundary is equal to the Bondi mass 
accretion rate {Mb = 4.6 x 10^^ g/s) for time durations of 
few X 10* yrs, amounting to few times the Bondi time. The 
duration for which Bondi accretion is reasonably reproduced 
lengthens with increasing outer radius of the computational 
volume. On adding an external galaxy bulge potential to the 
simulation, the mass inflow rate increases w.r.t. the Bondi 
rate, reaching a maximum value of I.GAIb- 

(3) The gas in our simulations accretes at the Bondi rate 
in the inner parts, between rin and rtr (the transition radius 
from net inflow inside rtr to net outflow outside it, i)3.1.4|) . 
but an outflow develops and dominates in the outer regions. 
We see that rtr = rout (i.e., the whole computational volume 
is accreting at the Bondi rate) only at very early simulation 
epochs. Later on, rtr reduces and starts to move inward. 

(4) More stringent boundary conditions can reduce the 
outflow at rout to a negligible fraction. Two of our tests were 
successful, where we change the velocity and pressure gra- 
dient (one in each case) of every particle having a net out- 
ward radial velocity {vr > 0) in a wide outer shell, between 
[0.4rout — rout]. With such a boundary condition, the mass 
inflow rate is nearly-steady (Min,ri„ — Mb) for durations up 
to ~ 25is. 



Next, we include radiative heating by a central X-ray 
corona and radiative cooling in our next series of simulations 
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(Table [2] 7run = 5/3), whose results are in the following 
(points 5 — 9): 

(5) The inner computational volume (r < 10 pc) is 
inflow-dominated for relatively-low values of the X-ray lumi- 
nosity {Lx /L-E&d = 5 X 10~^ — 5 X 10~^). In these simulations 
the particle velocity and density match the theoretical free- 
fall scalings well in the middle volume, with a deviation near 
the inner and outer radii, for the same reasons mentioned in 
point (1) before. The Mach number increases inward up to 
a certain r, but after that decreases at smaller r (because of 
a strong inner heating). 

(6) The runs with radiative heating and cooling take 
longer times to reach a steady-state, compared to the Bondi 
accretion runs. Some of the runs reach almost steady-state, 
when considering time intervals of few x 10* yrs. 

(7) In each run, some particles are in perfect radia- 
tive equilibrium. Some particles are above or below the 
T — ^ equilibrium curve, and are dominated by non-radiative 
heating and cooling, respectively. We identified three non- 
radiative mechanisms: heating near rin (at the highest-^, or 
lowest-r in a run) caused by adiabatic compression and nu- 
merical artificial viscosity, and adiabatic cooling near rout 
because of the outflow expansion at the outer boundary. 

(8) We find the occurrence of a heating of the inflowing 
gas in the inner volume, which is dominated by artificial vis- 
cosity. This is a powerful heating causing the temperature 
to rise to T ~ 10^ — 10** K, which is orders of magnitude 
above radiative equilibrium in some runs (Figure [T^ . Stan- 
dard AV prescriptions in SPH generates extra entropy at 
the expense of kinetic energy, therefore the velocity of the 
particles is smaller than the theoretical free-fall prediction 
at the same r in our simulations. The radius inside which 
the excessive AV heating occurs depends on the simulation 
parameters; we see it at r < 0.2 — 0.8 pc in different runs. 

(9) We detect signatures of radiative feedback between 
the characteristic X-ray luminosity of the central spherical 
corona (the radiation source, H2.4|) and the mass inflow rate 
at the inner boundary in our simulations. We see that the 
central infiow decreases at a varying rate as Lx is increased. 
High X-ray luminosities heat up the gas more, causing it 
to expand, increasing the outfiow at rout, and decreasing 
the inward flow at rin. With high enough values of Lx, the 
inflow at rin slows down, and an outflow develops in the 
inner volume. We see the outflow occurring for values of 
Lx/L^dd ^ 0.02, when = IQ''^^ g/cm^ 

We flnd non-spherical fragmentation and clumping of 
the gas in the runs with Lx/L-^dd = 0.01 and 0.02. This 
is likely due to a development of thermal instability which 
is triggered by small, non-spherical fluctuations in the flow 
properties inherent to SPH. We will present results of our 
detailed analysis of the clumping in a future publication. 
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APPENDIX A: SOLVING THE SPECIFIC 
INTERNAL ENERGY EQUATION 

We work out simple solutions of the energy equation inchxd- 
ing adiabatic and radiative terms, given radially-dependent 
velocity and density profiles. The hydrodynamic equation 
for the specific internal energy u is, 



(Al) 



Du 

'm \dt / p 

where, the first term in the RHS is adiabatic and the second 
is radiative. In terms of the temperature T — /i(7 — l)u/k, 
assuming steady-state (i.e., the time derivative is 0) dT/dt = 
0, for a radial velocity v = VrV, the energy equation can be 
rewritten as, 



dT 
dr 



(7-1) 



- — — (r'^v ) + -£. 

r2 fir ^ ^ ' h 



(A2) 



For a free-fall velocity (Ea. l2ip . the energy equation reduces 
to. 



dr 



= (7-1) 



ST 
27 



kv 



If 



(A3) 



We solve Eq. (|A3I) numerically, to find T as a function 
of r, using a fully explicit integration method with small 
enough Ar steps to prevent negative temperatures (Ar = 
0.002 - 0.0002 pc). 



dT 

r (r + Ar) = r (r) + (r) Ar. 

or 



(A4) 



An initial value of T = 2 x 10 K at r = rout — 200 pc is 
used during integration. We take the gas velocity as free- 
fall (Eq. [21) . and the free-fall density (Eq. [22)) assuming 
appropriate scaling from the particle snapshot plot. 
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Table 1. Simulations of Bondi Accretion 



Run 


"■out 


N ^ 


IC 








'end ° 


/IN 


/out 






h 














[Mq 


] 


r A -r 1 


[10 yrj 






[Mb 


(7ru 


n. Poo . Too 


M 


1 


5 


64^ 


Uniform ^ 




3.96 X 


10^ 


1.51 


3 


0.35 


0.65 






0.4 




2 


10 


64^ 


Uniform 




6.19 X 


lO" 


23.61 


7.2 


0.15 


0.85 






0.6 




3 


50 


128"^ 


Uniform 




7.73 X 


10* 


368.60 


20 


0.02 


0.9 






1 




4 


5 


64^ 


Bondi -i 




1.81 X 


10*5 


6.89 


2 


0.4 


0.55 






1 




5 


10 


64^ 


Bondi 




9.76 X 


lO" 


37.23 


8 


0.2 


0.8 






1 




6 


10 


128^ 


Bondi 




9.76 X 


10'' 


4.65 


8 


0.2 


0.8 






1 




7 


20 


128^ 


Bondi 




6.24 X 


10^ 


29.75 


8 


0.07 


0.9 






1 




7a k 


20 


128^ 


Bondi 




6.24 X 


10^ 


29.75 


80 


0.65 


0.05 






1 




7b > 


20 


128^ 


Bondi 




6.24 X 


10^ 


29.75 


100 


0.93 


0.01 






1.2 




8 


50 


128^ 


Bondi 




8.48 X 


lO'* 


404.35 


16 


0.02 


0.85 






1 




9 


20 


128^ 


PB. "init = 





6.24 X 


10^ 


29.75 


8 


0.06 


0.9 






1 




10 


20 


128^ 


Uniform 




4.95 X 


10^ 


23.60 


8 


0.06 


0.9 






1 




11 


20 


128^ 


Hernquist 




6.24 X 


10^ 


29.75 


7.2 


0.06 


0.85 






1.4 




12 " 


20 


128^ 


Bondi 




6.24 X 


10^ 


29.75 


8 


0.08 


0.8 






1.6 




" All th 


c runs 


have r^^ 


— 0.1 pc and 


7ru 


„ = 1.01. Init: 


ial conditions 


arc generated usin} 


? 7init 


= 1.01. p, 




io-i« g/ 




= 10^ 


K, and Ti„it = 


Too ■ Tlic corrcsp 


onding 


Bondi 


solution has 


Rb = 3.0 


pc, Rs — 


1.5 pc, 


and tB - 


= 7.9 X 10^ yr. 





^ N — Number of particles in the initial condition. 

^^tot.IC — Initial total gas mass within simulation volume. 

A/part — Particle mass. 
° ^cnd — Simulation end time. 

^ /in — Mass fraction accreted into r^^ by t — tend' 

S /out — Mass fraction moved outside r^^t t»y £ — ■ 

The steady-state value of mass inflow rate at the inner radius (if a steady-state is reached in a run), or the maximum 
value (if a steady-state is not reached by tend)' 

^ Uniform initial condition: Pi„it — Poo ■ nit ~ ^■ 

^ Bondi initial condition: Pinit(^) — PB ''^init — 

^ R,un with different outer boundary condition: Set Vj- — for particles which have Vj- > in between S — 20 pc. 

^ Run with different outer boundary condition: Set zero pressure gradient, {'VP/p)^ — 0, for all particles between 18 — 20 
pc and for those which have Vj- > in between S — 20 pc. 

Hernquist initial condition: Pinit(^) ~ P/j(^)i^init ~ 0- 

Run with a bulge potential of a Milky- Way type galaxy following the Hernquist profile fEq. llQl with iV-/j-,^]g(, — 3.4x10"'""'" Mq 
and a]3^igg — 700 pc. 
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Table 2. Simulations of Spherical Accretion with Radiative Heating and Cooling ^ 



Run 


rout 


AT 


MtatlC 




7init 


Too 


Rb 


Poo 


Tinit 


Lx 


tend 


/in 


/out 


A^iir.rjrr 


No. 


[pc] 




[Mq] 


IMq] 




[K] 


[pc] 


[g/cm3] 




[^Edd] 


[lO'S yr] 






[MB{7run,p'c^,Too)] 


13 


20 


128^ 


5.81 X 10^ 


0.277 


1.4 


10^ 


2.19 


10-21 


Too 


0.5 


1.0 


0.005 


0.99 


2 


14 


50 


128^ 


8.23 X 10'' 


3.92 


1.4 


lO" 


2.19 


10-21 


Too 


0.5 


2.9 


0.001 


0.99 


1.2 


15 


20 


128^ 


5.81 X 


2.77 X lO"'' 


1.4 


lO" 


2.19 


10-27 


Too 


0.5 


1.0 


0.005 


0.99 


2 


16 


20 


256^ 


5.81 X 10"^ 


3.46 X 10"^ 


1.4 


lO" 


2.19 


10-2" 


Too 


5 X 10-* 


1.9 


0.0052 


0.98 


2 


17 


20 


128^ 


5.81 X 10^ 


0.277 


1.4 


lO" 


2.19 


10-21 


Trad 1^ 


5 X 10-* 


2.9 


0.97 


0.03 


80 


18 


20 


128^ 


5.65 X 10^ 


0.269 


5/3 


10^ 


1.84 


10-21 


Trad 


5 X 10-* 


3.0 


0.97 


0.03 


80 


19 


20 


128^ 


1.47 X 10^ 


7.0 


5/3 


10^ 


183.9 


10-21 


Trad 


5 X 10-* 


1.5 


0.99 





2 


20 


200 


256^ 


1.33 X 10^ 


79.09 


5/3 


10^ 


183.9 


10-21 


Trad 


5 X 10-* 


6.5 


0.13 





6 


21 


200 


256"^ 


4 Q5 V 1 n'^ 


29.50 


5 /3 


10~ 


1.84 


10-21 


^ rad 


5 X 10 * 


8.7 


0.1 


0.012 


2000 


22 


200 


128^ 


1.33 X 10^ 


6.33 


5/3 


10^ 


183.9 


10-23 


Trad 


5 X 10"* 


70 


0.9 


0.1 


5 


23 


200 


256^ 


1.33 X 10^ 


0.791 


5/3 


10^ 


183.9 


10"23 


Trad 


5 X 10"* 


20 


0.3 


0.07 


4 


24 " 


200 


1.24 X 10^ 


9.77 X lO" 


0.791 


5/3 


10^ 


183.9 


10-23 


TRuir23 


5 X 10-^ 


19 


0.12 


0.01 


4 


25 


200 


1.24 X 10^ 


9.77 X 10** 


0.791 


5/3 


10^ 


183.9 


10"23 


TRuir23 


5 X 10"3 


21 


0.15 


0.02 


3.8 


26 


200 


1.24 X 10^ 


9.77 X 10'' 


0.791 


5/3 


10^ 


183.9 


10-23 


TRuir23 


1 X 10-2 


22 


0.2 


0.05 


3.4 


27 


200 


1.24 X 10^ 


9.77 X 10'' 


0.791 


5/3 


10^ 


183.9 


10-23 


TRuir23 


2 X 10~2 


25 


0.03 


0.1 


0.002 


28 


200 


1.24 X 10^ 


9.77 X lO" 


0.791 


5/3 


10^ 


183.9 


10"23 


TRuir23 


5 X 10-2 


50 


0.008 


0.99 


0.0001 



^ All the runs have r^^ — 0.1 pc and 'Yrun — 5/3. Initial conditions arc generated using Pinit(^) ~ and 'I'init C'") — ^s(^)- 

^rad ^ Initial temperatures set from the equilibrium T — § radiative heating-cooling function, i.e., the T solution obtained by solving for — in Eq. J13t . assuming 7 — 5/3. 
Runs 24 — 28 are started using the particle states from Run 23 at a time — 1.4 Myr as the initial condition, and changing Lx ii^ each case. 
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Table 3. Relevant Time Range for some Simulations 



Run 
No. 


fe/sl 


Is 
flO"* yrl 

L ^ J 1 


Tout 


Mtot,ic 


Mtot,ic/Ms " 
flO"* yrl 


N 


rout/(iV^) ^ 

LJr '-'J 


J. c 

f'Stcadv 

[10* yrl 


4 


4.6 X 10^^ 


0.79 


5 


1.81 X 10^ 


2.48 


64^ 


0.078 


0.6 


5 
(j 


4.6 X lO'-'^ 
l.G X 


0.79 
0.79 


10 

iO 


9.76 X 10^' 
9.76 X iO" 


1.3.38 
i:i.;i8 


64=^ 
128-'' 


0.156 
0.078 


1.6 
1.6 


7 


4.6 X 10^^ 


0.79 


20 


6.24 X 10^ 


85.52 


128^ 


0.156 


4 


8 


4.6 X 10^'^ 


0.79 


50 


8.48 X 10® 


1162.18 


128^ 


0.391 


8 



^ Mb = Bondi accretion rate for the simulation parameters: 7run = 1.01, poo = 10 g/cm^, Too = lO'^ K, 
Rb = 3.0 pc, Rs = 1.5 pc. 

^ tB = Bondi time. 

Mtot.ic / Mb = Time within which all mass would be accreted in, assuming accretion at the Bondi rate. 

rout/iN^^^) = Initial average inter-particle spacing, assuming uniform distribution of particles. It is a rough 
measure of simulation length resolution. Of course in the actual simulation the spatial resolution is given by the 
smoothing lengths of gas particles, which are dynamically determined at every timestep. 

isteady = Time for which simulation reproduces steady-state Bondi accretion. 
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r (pc) r (pc) 



Figure 1. Properties of particles in Run 7 as a function of radius, at simulation time t = 2tg = 1.6 X 10* yr. The quantities plotted 
from top-left: radial velocity, temperature, Mach number, density, acceleration, and smoothing length. The solid line is the median value, 
the dashed lines are the 95 percentile of the distribution lying below and above the median, and the dotted lines are the minimum and 
maximum values of the whole particle distribution at a given radius. The red curves in the top-left and middle-right panels show the 
Bondi solution for velocity and density. The horizontal straight lines in the top-left and middle- left panels indicate Vr = and Mach = 1. 



© 0000 RAS, MNRAS 000, 000-000 



18 P. Barai et al. 



/-476 





^j^j^^ 




-19 



-20 



-21 



-22 



Figure 2. Cross-section slice of gas density in the x 
with the velocity vector arrows. 



y plane through 2 = of Run 4 at time t = 0.25tg = 0.2 X 10'* yr, overplotted 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 
time (10" yr) 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 
time (10" yr) 
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time (10" yr) 
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o 
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/ 
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Figure 3. Mass inflow rate at the inner boundary as a function of time for the first eight runs in TablelTI Each panel has a different outer 
radius: rout = 5 (top-left), 10 (top-right), 20 (bottom-left), 50 pc (bottom-right), as the time coverage becomes longer. The top-right 
panel shows different particle numbers: N = 64'^ (Run 05) and 128^ (Run 06) for the Bondi IC. In addition, the top row and bottom-right 
panels show the results of the Bondi IC (Runs 04, 05, 06, 08), together with the uniform IC runs (Runs 01, 02, 03). The Bondi mass 
accretion rate (marked as the dash-dot-dot-dot horizontal line in each panel) is reproduced for a limited time duration. 
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2 3 4 
time (10" yr) 



Figure 4. Mass inflow rate at the inner radius vs. time for the last four runs in Table [T] The top panel shows results of different initial 
density profiles (see Table[T]and i|3.1.2l for details). The bottom panel overplots two runs with the same Bondi IC, but including a galaxy 
bulge potential in one case (Run 12), and not in the other (Run 07). 
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CO 
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Figure 5. Mass flux as a function of radius for four runs with different outer radius: rout = 5pc (top-left, Run 04), 10 pc (top-right, 
Run 06), 20 pc (bottom-left. Run 07), and 50 pc (bottom-right. Run 08). The inflow (solid line), outflow (dashed line), and net (dotted 
line) mass fluxes, as defined by Eq. I I20I I. are separately plotted. The absolute value of a mass fiux is shown whenever it is negative. 
is negative at all radii, and Mnct is negative in the inner parts where |A/in| > | A/out |. The simulation time of each run is noted in each 
panel. 
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Figure 6. Mass inflow rate at the inner boundary vs. time for the two runs in Table[T]with different outer boundary conditions, which 
have substantially reduced outflow (see ^3. 1.51 for details). Inflow is maintained for durations > 10 times larger than in Run07 in Fig.|3] 




Figure 7. Properties of particles in Run 23 (with heating and cooling) vs. radius at t = 1 Myr. The quantities plotted are: radial velocity, 
density, temperature and Mach number. The solid line is the median value, the dashed lines are the 95 percentile of the distribution 
lying below and above the median, and the dotted lines are the minimum and maximum values of the whole particle distribution at a 
given radius. The red curve in each panel shows the free-fall scaling of the corresponding quantity (discussed in il3. 2.111 . The blue curve 
in the first three panels depict the corresponding ZEUS results; dzeuS a-nd ^ff '^^'^ indistinguishable, pzEUS is at a constant offset from 
Pff because of different normalization. The green curve in the bottom- left panel indicate the free-fall temperature with only adiabatic 
processes. The horizontal straight lines in the top-left and bottom-right panels indicate Vr = and Mach = 1. 
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time (10'* yr) time (lO'^yr) 

Figure 8. Mass inflow rate at the inner boundary vs. time, for seven runs in Table[2]with Tinit = r^adi i-^-i the initial temperatures are 
in radiative equilibrium. Top-left panel shows that the accretion flow reaches almost steady-state with a larger outer radius. Top-right 
panel shows the effect of varying Too and poo . Bottom-left panel shows that different 7init has a negligible effect on the mass inflow. Note 
that the y-axis range in top and bottom panels are different. The details are discussed in ^3.2.21 



3 10' 



o 
to 




D5 10 



(D 
DC 



2 3 4 
time (10" yr) 

Figure 9. Mass inflow rate across the inner radius as a function of time for the flrst four runs in Table[2]with Ti^it = Too = 10^ K. The 
top panel shows the effect of varying outer radius (details in i|3. 2.21 1. and the bottom panel shows the results with different Lx ( ^3.2.41 1. 
The dash-dotted horizontal line marks the corresponding Bondi accretion rate. 
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1 10 100 1000 10000 0.1 1.0 10.0 

^ (erg cm / s) r (pc) 

Figure 10. Time evolution of particles in Run 18 (wliose mass inflow rate is shown in Figure |8j in the T — ^ and T — r planes. Each 
color represents the particles at a certain simulation time, which are labeled in the left panel in Myrs. Note that high ^ corresponds to 
small r. The left panel shows the median values of T and § at a given radius. The right panel shows the mean temperature and standard 
deviation around it plotted as vertical lines, indicating the particle distribution and scatter. 
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Figure 11. Temperature vs. photoionization parameter of particles (details in i|3.2.3l l. at a simulation time as labeled for each run, for 
the first four runs in Table [2] with Ti^it = lO'^ K (whose mass inflow rates are shown in Figure [SJ. The black line shows the radiative 
equilibrium T — ^ curve for the implemented heating-cooling function {C = 0, Eg. I13I I. The median values of T and J at a given radius 
are plotted as plus symbols, the maximum values as diamonds, and the minimum values as triangles. 
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Figure 12. Temperature vs. photoionization parameter of particles (discussed in ^3.2.311 . for nine runs in Table [2] with Tj^it = T^ad 
(whose mass inflow rates are shown in Figures [Sl and [T3)l . The median values of T and ^ at a given radius are plotted as the plus symbols. 
The radiative equilibrium T — § curve is overplotted as the red line in each panel. 
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Figure 13. Mass inflow rate at the inner boundary vs. time, for tiie last six runs in Table [2]with varying -Lx/iEdd between 5 X 10~^ 
and 5 x IQ-^. There is a significant drop in the mass inflow rate when Lx is increased from 0.01 to 0.02, which is associated with the 
development of an outflow. The details are discussed in i|3.2.4l 
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Figure 14. Radial temperature proflles of particles in three runs with different X-ray luminosity: X/x/i'Edd = 5 X 10~^ (Run 24), 
5 X lO"** (Run 23), 5 X 10"'^ (Run 25), at a simulation time t = 1.6 Myr. The solid line is the median value, the dashed lines are the 
95 percentile of the distribution lying below and above the median, and the dotted lines are the minimum and maximum values of the 
whole particle distribution at a given radius. The overplotted colored lines are the solutions of the steady-state internal energy equation, 
assuming free-fall scalings of the velocity and density: considering both adiabatic and radiative terms Tff ^r (EJ is plotted as the red 
curve, and considering only adiabatic terms Tfj ^ (Eg. I23II is plotted as the green curve. The details are discussed in i|3.2.5l 
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Figure 16. Properties of particles in Run 28 {Lx / L^dd = 0.05) vs. radius at i = 1.5 Myr, using a similar format as in Figure [7] The 
median values at few inner radial bins are plotted as plus symbols. 
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